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Figure 0: Ensemble averages of molecular dynamics simulations of laser-assisted particle removal. Images are line integral 
convolutions of averaged fluid vector field data from the clustered set of simulations where particle removal does not occur 
(left), from the entire data set of 100 simulations (middle), and from the clustered set of simulations where particle removal 
does occur (right). The particle is circular and rests upon a substrate, both shown as dark gray. Note that in the left image, 
where particle removal does not occur, the flow field at the bottom of the particle is directed around the particle, but this 
behavior is not as pronounced in the middle image, where all data sets have been averaged. Visible in the right image, beneath 
the particle, is the vertically directed flow field that is acting to lift the particle. 



Abstract 

Stochastic simulations of fluid flow deal with the evolution of un- 
certainties in initial and boundary conditions, parameters, and phys- 
ical models. These uncertainties may lead to qualitatively different 
behaviors of the flow, raising the question of how to visualize ran- 
dom flow variables in a meaningful way. Our approach is to cluster 
the simulation data by qualitative phenomenon and then to com- 
pute average flow quantities within each cluster, a technique we 
call "clustered ensemble averaging." This tactic extends the basic 
visualization strategy of extracting features from within a computa- 
tional domain defined within time and space. The abstract features 
we seek are defined across the space of all possible simulations: 
each feature (cluster) is the subspace of simulations that share a ba- 
sic behavior, We illustrate this technique on data from molecular 
dynamics simulations of laser-assisted particle removal, where the 
explosive evaporation of a laser-heated fluid acts on minute particle 
contaminants on a substrate, sometimes removing them and some- 
times failing to do so. We ran the simulation a hundred times, com- 
prising a hundred samples of the simulation space, and clustered 
each according to the behavior it exhibits. The images that result 
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from clustered ensemble averaging reveal characteristics common 
to each cluster that are hidden when the average is calculated over 
all outcomes. 
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1 Introduction 

One of the main strategies of visualization is to display salient fea- 
tures in a dataset so they can be interpreted by the scientist, engi- 
neer, or physician who is analyzing the data. For three-dimensional 
(3D) scalar fields, feature-extraction techniques have been devel- 
oped to display boundaries of tissue in medical data [7, 18]. For 
3D vector fields, feature-extraction techniques have been developed 
to display shock fronts [15] and vortex cores [3, 12, 16]. Some 



of these techniques have been extended to extract features from 
4D time-varying datasets, where the feature extends across both 
space and time [5, 19, 23, 25]. Additional techniques capture fea- 
tures that exist in 4D scale space, where cross-sections at differ- 
ent values of the scale axis correspond to different resolutions of 
the data [10, 14, 24]. These examples illustrate feature-detection 
techniques that have been developed to locate and display features 
in two or three spatial dimensions, and in multidimensional do- 
mains of space-time or scale space. Other visualization techniques 
merely approximate the geometric appearance of a feature, match- 
ing its shape to a best-fit model such as an ellipsoid [17]. This 
approach can be understood as finding an average approximating 
shape among all features in a certain equivalence class. 

We extend the idea of visualization via feature-detection to a new 
level of abstraction by considering a feature that exists across not 
only in time and space, but also across a set of variations. These 
variations arise from multiple runs of a stochastic (random) simu- 
lation starting from different initial conditions. 

The three-body problem is a simple example that illustrates how 
diverse outcomes can result from minute variations of initial condi- 
tions: in the plane, three particles of equal mass are assigned posi- 
tions pj,p 2 ,p 3 and velocities Vj^jVj at time / = 0. Each particle 
has four independent parameters (two for position and two for ve- 
locity), so the planar three-body problem is characterized by a 12- 
dimensional parameter space. Although most solutions of the pla- 
nar three-body problem appear random and are of minimal interest, 
occasionally a periodic solution is obtained. Finding parameters 
that produce the the qualitative behavior of periodicity is still an 
area of active research [8]. In comparison, a computational molec- 
ular dynamics simulation of a thin layer of fluid may have more 
than a thousand particles. Although the 4000-dimensional parame- 
ter space of initial conditions is very large, the number of qualita- 
tively distinct outcomes may be small. 

Stochastic processes may be intuitively defined as systems which 
evolve probabilistically in time. Mathematically described, a 
stochastic process is a time-dependent random variable or field that 
evolves in time. In the present context, these processes are obtained 
as solutions to evolution equations with random initial and bound- 
ary conditions, and/or random physical and geometric parameters. 
The stochastic simulations of a flow configuration (created by some 
set of external conditions) are similar in some sense but differ from 
one another in detail. In order to analyze the underlying phenom- 
ena in such flows, one needs to define the mean value of the flow 
properties. A meaningful mean value is defined by the statistical 
average of the ensemble of all similar flows. This permits decom- 
position of a field variable into a mean and a fluctuation (deviation). 
The mean field then allows one to extract the "deterministic" orga- 
nized structures underlying the random flow. The visualization of 
probabilistic evolution of fluid flows can be considered a form of 
"uncertainty visualization" [9, 11]. 

We present a technique for displaying the mean behaviors of a 
flow that exhibits two qualitatively distinct outcomes. Multiple runs 
of the simulation yield datasets that are assigned to one of two clus- 
ters. Ensemble averaging within one cluster reveals behavior that 
is not evident from taking an overall average. We demonstrate the 
technique on an example: a simulation of laser-assisted particle re- 
moval. 

The remainder of the paper is organized as follows. Section 2 
provides background into laser-assisted particle removal. Section 3 
describes the computational simulation model. Section 4 describes 
how the particle data from the simulation is converted into contin- 
uous density fields and how these datasets are clustered. Section 
5 describes how the particle data from the simulation is converted 
into continuous velocity data. Section 6 summarizes the results. 
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Figure 1 : Laser assisted particle removal (LAPR). Left panel shows 
circular particle on substrate surface. Vaporized energy transfer 
medium (ETM) is directed toward the substrate surface. The ETM 
condenses on the surface as shown in the middle panel. Laser en- 
ergy is then directed toward the particle and liquid ETM, causing 
explosive boiling of the ETM and removal of the particle, as shown 
in the right panel. 



2 Laser-Assisted Particle Removal 

Manufacturing processes for hi-tech devices in microelectronics 
or data storage continue to demand higher resolution, smaller 
linewidths, miniaturization, and vanishing mechanical clearances. 
Since it is not currently possible to totally eliminate all contamina- 
tion sources during industrial processing, there is an ever-increasing 
need to develop methods for removing minute particle contami- 
nation from critical surfaces. Currently, the most widely used in- 
dustrial cleaning techniques use liquid chemicals, which are them- 
selves sources of contamination. Removing these particle contami- 
nants is difficult since a submicron particle can adhere to a surface 
with a force that is over a million times its weight [13]. A cleaning 
technique for submicron particles should be effective for particle 
removal, not cause surface damage, or add any extra contamination 
to the surface. Additionally, we want the process to be efficient, 
simple, fast and chlorofluorocarbon-free. 

Laser assisted particle removal (LAPR) is a technique, which 
meets the criteria listed above, for removing small particles (as 
small as 50 nm, though currently used principally for particles 
~ \fim in diameter) from a substrate. In this process, a small 
amount of energy transfer medium (ETM), usually a mixture of 
water and alcohol, is condensed onto the substrate, and then laser 
energy is directed at the substrate surface, as shown in Figure 1. 
This energy is absorbed by the ETM fluid, and perhaps also by the 
substrate and/or the particle, and results in removal of the ETM 
and sometimes in the removal of the particle. Particle removal is 
thought to be caused by a combination of three mechanisms: ex- 
plosive boiling of the liquid ETM, thermal expansion of the particle 
("hopping" effect), and thermal expansion of the substrate ("tram- 
poline" effect) [22]. 

A single optimized laser pulse has been shown to effectively re- 
move 90% of particles with sizes on the order of 1 Jim [2]. Re- 
peating the process allows removal efficiency of nearly 100% to be 
achieved, typically within five to eight repetitions [2]. Two major 
factors influence the particle removal efficiency of LAPR: (1) the 
adhesion forces holding the particle to the substrate surface, and 
(2) the laser-induced particle removal forces. It is well known that 
the cleaning efficiency increases with increasing laser fluence, but 
at very high laser fluences substrate surfaces are easily damaged 
by laser irradiation. Thus, determining the optimal laser cleaning 
conditions and clearly understanding the interaction mechanisms 
between particle and substrate surface are the goals of modeling 
the LAPR process. 



3 Simulation model 

A two-dimensional molecular model was constructed for simulat- 
ing LAPR. The model consists of a substrate, a particle, and the 



Table 1 : Lennard-Jones potential parameters for interactions be- 
tween the ETM, particle and substrate. 
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ETM fluid. To reduce computational requirements, a very simple 
molecular model was applied. All interactions between molecules 
in the ETM fluid, substrate, and particle were modeled using the 
Lennard-Jones 12-6 potential [1]: 



Vij) i r u 



(1) 



where r t j is the distance between molecules i and j, G is a measure 
of the molecule's diameter (the distance where the energy is zero), 
and -e is the minimum value of the energy. The simulations were 
run using reduced (nondimensional) units [1], with the mass of all 
molecules equal. (For convenience, in addition to reduced units, 
we will present dimensional units assuming argon molecules; For 
argon, a = 0.43 nm, and £ - 1.67 x 10~ 21 J.) The Lennard-Jones 
potential parameters for the interactions between the ETM fluid, 
particle, and substrate are given in Table 1 . 

The substrate is composed of five layers of molecules arranged in 
a hexagonal lattice, which is the two-dimensional crystal structure 
of the Lennard-Jones solid. During the simulation, the molecules of 
the substrate interact (via the potential function) with nearby fluid 
and particle molecules, but the substrate molecules are held in fixed 
positions. This is partly for additional computational cost savings, 
but mainly because thermal expansion of the substrate in the peri- 
odic geometry would likely cause buckling, which would substan- 
tially complicate our analysis. The rigid substrate approximation 
avoids this, although it does mean that the "trampoline" mechanism 
for particle removal is not present in these simulations. 

The particle is circular, with a diameter of 19<J (6.46 nm), and 
is also composed of molecules in a hexagonal arrangement. The 
lowest layer of molecules that compose the particle is in contact 
with the molecules in the uppermost layer of the substrate. 

Simulations were conducted on a two-dimensional computa- 
tional domain of size* = 202a (68.68 nm) byy = 240a (81.6 nm). 
The domain is periodic in jc and molecules that reach the maximum 
y-value are reflected back into the domain. This reflection is nec- 
essary to prevent the entire fluid from evaporating through the top 
of the domain. Since the important part of the simulation occurs in 
the lower half of the domain, the reflection condition has no effect 
on the results that follow. 

After placing a particle on the substrate within a computational 
cell, fluid is added to the cell by performing a Grand Canonical 
Monte Carlo (GCMC) simulation. GCMC is a stochastic technique 
in which the number of molecules is allowed to vary according to 
a specified chemical potential [1]. This is accomplished through 
"moves" in which molecules are either displaced, created at ran- 
dom positions, or destroyed. When completed, as determined by 
stabilization of the total number of molecules, a liquid (ETM) film 
of varying thickness overlies the substrate and the particle, as shown 
in the first column of Figure 2. 

Because particle removal is a discrete event, in order to deter- 
mine removal efficiency it is necessary to simulate the removal pro- 
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Figure 2: Snapshots from LAPR simulations. The substrate is 
shown as a dark gray band at the bottom of each panel, the ETM 
fluid molecules are black, and the particle, which starts in contact 
with the substrate, is circular and gray. Each row represents a sim- 
ulation that began with slightly different initial conditions. Each 
column is time advanced by 4,000 iterations (0.044 ns). Substrate 
temperature is 1. 7 (205. 7°K). ETM fluid thickness is 50c (1 7 nm). 
In the right column, we see that particle ejection occurs in half of 
the simulations. 



cess many times, starting from various initial conditions. Addi- 
tional configurations were generated by allowing the GCMC simu- 
lation to run for an additional short time (100 cycles). Initial con- 
figurations for an ETM thickness of 50a (1 7 nm) are shown in the 
leftmost column of Figure 2. 

In the molecular dynamics simulations, a fifth-order Gear 
predictor-corrector algorithm was used to integrate the equations 
of motion, with a time step of 0.005 (5.5 x 10" " s). After 250 time 
steps, the substrate temperature was instantaneously increased from 
its base value of 0.4 (48.4° K) to a value between 0.5 (60.5°K) and 
5.0 (605° K). Increasing the substrate temperature in this manner 
approximates the effect of laser heating in the case where the laser 
energy is absorbed by the substrate, rather than by the fluid. For 
more details on the LAPR simulations, see [2 1 ]. 

Each simulation was run for 30,000 time steps (representing 
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Figure 3 : Snapshots of LAPR simulations. For all simulations, the 
substrate temperature, T, is 1.7 (205. TK) and the elapsed time, t, 
is 30,000 (0.33 ns). Each column shows a different ETM thickness, 
W, with W increasing to the right, and taking on the values 6a 
(2.04 nm), 10c (3.4 nm), 25o (8.5 nm), 50o (17 nm), and 70a 
(23.8 nm). Each row shows a simulation that began with slightly 
different initial conditions, X. 



0.33 ns of physical time), which, except in the case of tempera- 
ture equal to 0.5 (60.5°K), was enough to observe the complete 
ejection of the ETM fluid layer, and, in some cases, particle re- 
moval. Typical results are shown in Figure 2. In this figure, each 
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Figure 4: Averages of snapshots of LAPR simulations. Elapsed time 
t is 30,000 (0.33 ns). ETM thickness W increases to the right, and 
substrate temperature T increases upward, from 1.5 (181.5° K) to 
2.0 (242° K) by 0.1 (12.1°K). Note that the combination of high tem- 
perature and medium width (upper row, middle column) results, on 
average, in particle removal. 



column is time advanced by 4,000 iterations (0.044 ns), and each 
row began with slightly different initial conditions (i.e., number of 
ETM molecules, position of ETM molecules, and velocity of ETM 
molecules). The substrate temperature for each of these simulations 
was 1.7 (205. 7° K). By examining the rightmost column, we see that 
of the ten simulations, five resulted in particle removal, and five did 
not. Since the initial conditions for the set of ten simulations varied 
only slightly, and the substrate temperature was identical for each 
simulation, why was the particle removed in some simulations and 
not in others? 



4 Visualization 

Before describing the process of visualizing the fluid flow to an- 
swer this question of why the particle is sometimes removed, we 
begin by defining variables. Let T be the substrate temperature, t 
be the elapsed time, W be the width of the ETM fluid layer, and 
X represent the simulation number (or the initial conditions of the 
simulation). Then, Figure 2 shows / increasing to the right and X 
increasing upward, while T ~ 1 .7, and W = 50o\ 

We can also fix the elapsed time, t, and the temperature, F, and 
vary the initial conditions, X 7 and the ETM thickness, W, as shown 
in Figure 3. This figure consists of snapshots from LAPR simula- 
tions at t = 30,000 (0.33 ns), and T = 1.7 (205.7°K), where in this 
case, the ETM fluid thickness, W, varies from very thin (6a or 2.04 
nm), to very thick (70a or 23.8 nm). Again, notice the range of 
behaviors that occur in LAPR: sometimes the particle lifts with the 
ETM fluid, and sometimes it doesn't, even as the thickness of the 
fluid layer varies. Looking at each snapshot in this manner, it is 
difficult to extract meaningful data from the simulations since the 
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Figure 5: Ensemble averages of LAPR simulation snapshots for 
simulations where there is not particle removal (upper) and for sim- 
ulations where there is particle removal (lower). Elapsed time t is 
30,000 (0.33 ns). ETM thickness W increases to the right, and sub- 
strate temperature T increases upward, from 1.5 (181.5° K) to 2.0 
(242°K). 



Figure 6: Averages of 10 (left) and 100 (right) simulations at tem- 
perature T = 1.6 (193.6° K), elapsed time t = 30,000 (0.33 ns), and 
ETM thickness W = 50a (17 nm). 



1 



Figure 7: Ensemble averages (over time t and simulation X) of par- 
ticle trajectories for the case of no particle removal (left) and par- 
ticle removal (right). T - 1.6 (193.6° K), and W = 50g (17 nm). 
(The contrast was adjusted to enhance the image.) 



dimension is so large. 

We get a better picture of the average response for a given T and 
W by averaging over different initial conditions; that is, we integrate 
out the variable X. This gives rise to Figure 4, which shows fluid 
thickness, W, increasing to the right, and temperature, 7\ increasing 
upward from 1.5(181.5°K) to 2.0 (242°K),by 0.1 (12.1°K). In par- 
ticular, the third row from the bottom shows the results of averaging 
out all the X's from Figure 3. Notice there is a bimodal appearance 
of the particle, especially in the middle columns, where you can 
see an average position of the particle in contact with the substrate 
and another average position between the substrate and the average 
lifted fluid layer. Each of these two modes corresponds to a subset 
(called an ensemble) of the simulations shown in Figure 3. Rather 
than averaging all of the simulations, we can instead average only 
those simulations in the first ensemble, where the particle is not 
removed, as shown in the upper panel of Figure 5, and then aver- 
age the simulations in the second ensemble, where the particle is 
removed, as shown in the lower panel of Figure 5, The missing en- 
tries (empty rectangles) indicate that there is no data. For example, 
at the bottom right of the lift ensemble (lower panel) we see that 
particle removal did not occur in any of the simulations. 

For the simulations presented thus far, there were 10 initial con- 
ditions, X, for each ETM thickness, W, and temperature, T. For 
comparison, we also conducted 90 additional simulations for the 
case T = 1.6 (193.6°K), and W = 50a (17 nm), making a total of 
100 simulations for that case. Figure 6 shows the average of the 
initial 10 simulations (left panel) and the average of all 100 simu- 
lations (right panel) at / = 30,000 (0.33 ns). With more simulations 
to average, the variance is reduced and the resulting average image 
looks much smoother, but at the expense of more compute time. 
We note that each simulation took about 14 minutes on a desktop 
computer with a 1.2 GHz processor and 256 MB of RAM. 

Although 100 simulations provides a smoother average, the tra- 
jectory of the particle is difficult to perceive. By calculating the en- 



semble averages of only the particle molecules (averaged over both 
time, t y and initial conditions, X), the average particle trajectories 
are easily seen, as in Figure 7. On the left, is the ensemble average 
trajectory of the particle for the 30 simulations where particle re- 
moval did not occur. In this case, it can be seen that sometimes the 
particle rolls to the left or right, but the most likely (darkest) case is 
no motion at all. On the right is the ensemble average trajectory of 
the particle for the 70 simulations where particle removal did occur. 
In this case, it can be seen that the most likely (darkest) case is that 
the particle lifts off the substrate perpendicular to the substrate, al- 
though it is also possible that the particle will lift-off the substrate 
at a slightly different angle. 

A faster way to produce a smoothed image of the molecular den- 
sity (as shown in Figure 6) is illustrated in the upper panels of Fig- 
ure 8. The upper left panel shows the molecules in space with a grid 
laid over the domain. Each molecule is assigned a density function 



y;.(x) = 5(x-x | .), 



(2) 



where p t - is the position of / and S is the Dirac delta distribu- 
tion. We convolve / =51/) with a filter function h(x) to blur the 
molecular densities across the image, as shown in the upper mid- 
dle panel of Figure 8, We choose a filter function with compact 
support, namely the cubic approximation to the Gaussian given by 
h(x) = 2JC 3 - 3JC 2 + 1 , where x is scaled so that x = 1 at the filter 
radius, which is given in spatial coordinates. We can think of the 
resulting image, /(x), as a scalar function from R 2 -¥ R, corre- 
sponding to the molecular density at a point in space, where / is 
given by 

£,./(x)ft(x- P| .) 

' W 2,A(*-P f ) ' C) 
and is illustrated in the upper right panel of Figure 8. We normalize 
(divide by the sum of h) because the superposition of filter functions 
may sum to something other than one. 

The lower panel of Figure 8 shows the result of this convolution 
on a given simulation, producing a sequence of images with dif- 
ferent filter radii. The radii used, from left to right, are la (0.34 
nm), 8a (2.72 nm), and 32a (10.88 nm). In each image the color 
black corresponds to low density (0 molecules per unit area), while 
white corresponds to high density (typically about 60 molecules per 
unit area). Thus, we have transformed a collection of about 8,000 
molecules into a continuous scalar-valued density function. 

The upper panel of Figure 9 shows smoothed density images 
for a time sequence of simulations where T - 1.5 (181.5°K), and 
W = 25 a (8.5 nm). If we treat the time direction as a third spatial 
coordinate, these density functions can be stacked up to produce 
a 3D scalar field /(*,>>,/) as shown in the lower left panel of Fig- 
ure 9. An isosurface of this scalar field with density value of about 
30 is shown in the lower right panel of Figure 9, where we see 
the boundary of the subvolume which contains the bulk of the fluid 
molecules. This figure also shows the trajectory of the particle (as 
a dark cylinder). The bottom of the trajectory corresponds to the 
beginning of the simulation where the particle is on the substrate. 
At the top of the trajectory, where t = 20, 000 (0.22 ns), the particle 
has moved away from the bottom of the domain, having been lifted 
by the ETM fluid which has moved past the particle. 

Figure 10 shows the two lifting behaviors more clearly. With in- 
creasing time, the isosurface containing the volume of fluid moves 
away from y = 0. As in Figure 9, the time direction is aligned ver- 
tically, and the y direction of the image is aligned in the horizontal 
right direction. The left panel shows a run where the particle lifts 
and is removed, while the right panel shows a run where the particle 
does not lift. In both panels, notice that the ETM fluid closest to the 
substrate (on the left in the figures) begins to move before the fluid 
farthest from the substrate (on the right in the figures). This differ- 
ence in velocities is visible as a kink in the isosurface near t = 0 on 
the right side of both figures. 
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Figure 8; Smoothed density images. Upper panels illustrate how 
smoothed density images are computed: left panel shows molecules 
in space with an overlaid grid, middle panel shows molecules each 
with its density function centered at the molecule 's position, and 
right panel shows the resulting density image, lower images are 
examples of smoothed density images with varying filter sizes. The 
filter radii used, from left to right, are la (0.34 nm), 8o (2. 72 nm) 
andS2o (10.88 nm). 




Figure 9: Smoothed density images are used to make an isosurface. 
Upper panel: time sequence of smoothed density images. Lower 
left panel: Smoothed density slices stacked in time. Lower right 
panel: Isosurface of density, with particle trajectory shown as dark 
cylinder. 




Figure 10: Left: Close-up of density isosurface (right panel of pre- 
vious figure) for a simulation where particle removal occurs. Right: 
Close-up of density isosurface for a simulation where particle re- 
moval does not occur. Time, t, increases in the vertical direction. 
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Figure 1 1 : Continuous velocity fields. Upper panels show how a 
continuous velocity field is constructed from molecular velocities. 
Upper left panel: molecules in space with a grid overlaid. Up- 
per middle panel: each molecule is shown with the filter function 
centered at the molecules position. Upper right panel: after sum- 
ming the product of the filter function with the velocity, and nor- 
malizing, the result is a continuous vector field on the grid. Lower 
panels show continuous vector fields using line integral convolu- 
tion. Images of continuous vector fields are for the case of T - 1.7 
(205. 7°K), W = 50a (1 7 nm), t - 30,000 (0.33 ns), andX = /, with 
different filters of varying widths. The widths of the filters used, 
from left to right, are la (0.34 nm), 8a (2. 72 nm), and 32a (10.88 
nm). 



5 Vector fields and ensembles 

The previous section showed ensemble averages that distinguish 
qualitative behavior of the particle. Figures ?? explicitly show ex- 
amples of particle trajectories in the two regimes (particle removal 
and no-removal), and of the ETM layer's trajectories in these two 
regimes. It is clear that two distinct behaviors occur; the question is, 
what property of the ETM causes the difference? The ETM always 
lifts whether the particle does or not, and at the crucial early stage 
of the simulation, the images of the ETM density are practically in- 
distinguishable. We were therefore motivated to look at clustered 
ensemble averages of the ETM velocity for evidence of any subtle 
difference between the two regimes. This required us to produce a 
vector field from discrete particles in motion. 

Just as we convolved molecular densities with a filter function 
to produce a continuous scalar density field, we can also convolve 
the molecular velocities with a filter function to produce a contin- 
uous velocity field. The upper panel of Figure 1 1 illustrates this 
process. The upper left panel shows several molecules each with a 
given velocity vector. A grid is laid over the domain. The upper 
middle panel of the figure shows the filter function centered at each 
of the molecular positions p t . The vector valued velocity at each 
grid location is found by summing the product of the filter function 
and the velocity for all molecules: 



(4) 



Again, we normalize by dividing by the sum of h since the sum 
of filter functions applied at a grid point may be greater than one. 
As before, we used the cubic interpolating polynomial h{x) = 
Ix 3 - 3X 2 + 1, for the filter function, where jc is scaled so that x 
= 1 at the filter radius, which is given in spatial coordinates. The 



final velocities on the grid are shown in the upper right panel of 
Figure 1 1 . 

The lower panels of Figure 1 1 show the continuous vector fields 
of a given simulation for three different filter radii using line inte- 
gral convolutions (LIC) [6] to display the vector fields. The LIC 
images show streamlines of the flow The solid black regions corre- 
spond to missing data where no molecule was close enough for its 
velocity to be blurred at that point. On the left, using a small radius 
of la (0.34 nm), the sparseness of the data is apparent. The mid- 
dle panel shows that by increasing the radius of the filter function 
to 8a (2.72 nm), the missing data can be filled in, leaving missing 
data mostly at the bottom of the image. In the right panel, we have 
further increased the filter radius to 32a (10.88 nm). We see that 
although it is possible to fill in the missing data making a dense 
vector field, much of this vector field corresponds to sparse data. 

Although the images in the lower panels of Figure 1 1 were cre- 
ated using only one simulation, the same process can be applied to 
many simulations creating an image of the average velocity over 
many runs (namely, the ensemble-averaged velocity). We see the 
results of such an average in the teaser image (Figure 0; at the be- 
ginning of the paper). For the center panel of Figure 0, we averaged 
100 simulations at T = 1.6 (193.6°K), t = 30,000 (0.33 ns), and W 
- 50a (17 nm). For the left image, we averaged the ensemble of 
simulations for the cluster where particle removal did not occur (30 
cases), and for the right image, we averaged the ensemble of simu- 
lations for the cluster where particle removal did occur (70 cases). 
In each of these images, we used a filter radius of 8a (2.72 nm), 
which corresponds to the lower middle panel of Figure 1 1. In each 
of these images, the particle and substrate are shown as gray re- 
gions to allow these locations to be seen against the vector field. 
In the left panel (the cluster that does not lift) the velocity field at 
the base of the particle is clearly directed in an outward direction 
forming what looks like a wake around the particle, while in the 
middle panel (all images averaged together) this behavior is not as 
pronounced. The center panel more closely matches the right panel 
where particle removal has occurred. In the right panel of Figure 0, 
in the space between the particle and the substrate, the velocity field 
beneath the particle can be seen acting on the particle in an upward 
direction, suggestive of the lifting force that successfully dislodges 
the particle from the substrate. 

6 Verification and summary 

Our research group includes two persons with prior involvement in 
simulation of laser-assisted particle removal (Smith and Hussaini). 
Although it seemed evident to us that the clustered ensemble aver- 
ages of the energy transfer medium (ETM) displayed different qual- 
itative behaviors in particle removal versus no-removal, we may 
have been biased due to our own involvement in developing this 
visualization tool. To validate the usefulness of clustered ensemble 
averaging as a visualization technique, we recruited an expert to an- 
alyze the resulting images. The expert, a chemist, has a laboratory 
to collect experimental data from actual LAPR trials, and so was 
familiar with the underlying science. He was not involved in the 
development of the visualization technique, and had not seen the 
images resulting from clustered ensemble averaging. 

The protocol for our simple user-study was very informal: we 
conducted a telephone interview with the expert. We first asked 
whether any observable difference was expected in the ETM fluid 
for the particle removal versus the no-removal case. We then 
pointed him to the images, which we had put on the World Wide 
Web, and answered any questions about how the images were pro- 
duced. We then asked for his interpretation of the images. 

The expert expected there to be observable differences in the 
ETM fluid for the particle removal versus the no-removal case, but 
was unable to articulate what that difference would be. He reported 



seeing clear differences in the clustered ensemble averages of the 
ETM velocity (Figure 0, left and right panels). He remarked that 
in the case of particle removal, "you can plainly see the laminar 
flow [vertical straight lines]." For the no-removal case he said, "the 
flow lines around the particle are like those you would see around a 
ball in a wind tunnel," referring to the appearance of the wake-like 
velocity structure around the particle. 

Clustered ensemble averaging reveals differences between be- 
haviors of the particle removal and no-removal regimes that an ex- 
pert chemist has confirmed. The visualization technique gave him 
new insight and raised new questions about the behavior of fluid. 
This indicates that clustered ensemble averaging can play an im- 
portant role in coupling stochastic simulation with experimental 
science by providing a visual tool for the analysis of qualitative 
behaviors. 
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